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SUMMARY 

The three-dimensional Navier-Stokes solution code using the LU-ADI factorization 
algorithm was employed to simulate the workshop test cases of transonic flow past a 
wing model in a wind tunnel and in free air. The effects of the tunnel walls are well 
demonstrated by the present simulations.. An Amdahl 1200 supercomputer having 128 
Mbytes main memory was used for these computations. 

INTRODUCTION 

Computational Fluid Dynamics (CFD) is achieving a status comparable to classical 
methods of laboratory experiment and theoretical analysis. The foremost tools of CFD 
researchers are supercomputers. With the available huge memory of these machines, it 
is now possible to solve the compressible Navier-Stokes equations using a large number 
of grid points. 

Since the use of a supercomputer is one of the important factors for practical Navier- 
Stokes simulations, the solution algorithm which is used should have good vectorization 
capability in addition to efficiency and robustness. Efforts have been made by the 
authors and many scientists to develop such algorithms, and a number of results using 
supercomputers have been reported recently [1-3]. The present authors also have been 
engaged in developing an efficient and robust three-dimensional Navier-Stokes code. 
The LU-ADI algorithm used in the code was originally developed by Obayashi and 
Kuwahara [4] for the two-dimensional problem. Extension to three dimensions and 
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application to relatively simple geometries were reported in reference 5. A series of 
studies [6-10] for practical applications of the code to transonic flow over transport 
aircraft configurations have been done in parallel with improvements in stability and 
efficiency. 

In the present paper, computations for the workshop test problems were carried out 
by using the three-dimensional Reynolds-averaged “thin-layer” Navier-Stokes equa- 
tions. The turbulent eddy viscosity was computed by using the Baldwin-Lomax model 
modified for two thin layers [11]; that is, for the wing surface and wind-tunnel walls. 

The LIJ-ADI factorization algorithm was used to obtain the steady-state solution. 
Computations were carried out for 0 and 2° angle-of-attack cases for the wing in the 
wind tunnel, and for 0, 2, 5 and 8° angle-of-attack cases for an isolated wing with the 
free-stream Mach number and Reynolds number corresponding to the experiment [12]. 
An Amdahl 1200 supercomputer was used for the computations. For the wind-tunnel 
simulation, the CPU time was 10.0 n sec/step/grid point, and the required memory 
was 92 Mbytes for about 0.8 million grid points. 

GOVERNING EQUATIONS AND NUMERICAL ALGORITHM 


Compressible Navier-Stokes Equations 

The partial differential equations governing the three-dimensional compressible vis- 
cous flow of an ideal gas are the compressible Navier-Stokes equations. Written here 
in the conservation-law form and in the generalized coordinates, they are given by, 

d T Q + d^E + d v F + d^G = Re~ 1 (d v S i + d f 5 2 ). (1) 

For high Reynolds-number flows, the viscous effects are confined to a thin layer near 
wall boundaries and are dominated by the viscous terms associated with the strain 
rates normal to the surface. The terms associated with the strain rates along the 
surface are comparatively small and negligible. Also, since the computational grid is 
highly concentrated near the surface and, as a result, the mesh aspect ratio becomes 
very large, inclusion of these terms would not change the solution. Thus, the thin-layer 
approximation has been introduced in equation (1) as is typically done in most Navier- 
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Stokes computations. In addition, the thin-layer approximation has been extended to 
two directions in equation (l) since viscous layers must be considered on both the wing 
and the wing root at a side wail. The cross-derivative terms are neglected and only 
and S 2 ? terms are retained. 

The pressure, density, and velocity components are related to the internal energy for 
an ideal gas as follows, 


p = {7 - l){e - p{u 2 + v 2 + w 2 )/ 2 ). (2) 

In the following computations, the viscosity coefficient in equation (l) is computed 
as the sum of (x + fx t . The turbulent eddy viscosity, fi t , is computed by using the two- 
layer algebraic eddy viscosity Baldwin-Lomax model with a modified length scale. The 
length scales are difficult to define near the junction of a fuselage and a wing. In the 
present study, the evaluation of length scales proposed by Hung and Buning for their 
blunt-fin study [11] is adopted. 

The metrics are evaluated using second-order central-difference formulas for interior- 
points, and three-point, one-sided formulas at the boundaries. The steady-state solu- 
tion of equation (1) is obtained by time integration in a time-asymptotic fashion. 

LU-ADI Factorization Algorithm 

The time-marching method used here is the LU-ADI factorization proposed by the 
present authors [8]. Implicit time-integration methods with a delta form are widely 
used for solving steady-state problems since the steady-state solution is independent of 
the left-hand side operators. The most commonly used is the ADI factorization Beam 
and Warming method used in reference 13. The present LU-ADI factorization method 
is a compromise of ADI and LU factorization [14]. Each ADI operator is decomposed 
into the product of lower and upper bidiagonal matrices by using a flux vector-splitting 
technique [15] and a diagonally dominant factorization [16]. 

It should be noted that the right-hand side of the equation remains in the form of 
central differencing and numerical smoothing similar to the Beam- Warming method. 
For the convective terms on the right-hand side, fourth-order differencing is used except 
near the boundaries. Maintenance of the free stream is achieved by subtracting the 
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free-stream fluxes from the governing equations. As a result, the present scheme is first 
order in time and second order in space. 

The present method applied to equation (1) is written as, 

T^aDaU^T^LbDbUbT^T.LcDcUcT- 1 ^ 

= + Sr,F n + 6;G n - Re-'(6%S J* + 6 2 S 2 n )} 

+(D e \ i + D t \ rl + D e \ i )Q n (3) 

where h is the time step, S is a central finite-difference operator, and D e is the explicit 
nonlinear artificial dissipation operator. On the left-hand side, the original ADI oper- 
ator of the Beam- Warming method is rewritten using the diagonal form [17] and the 
flux vector splitting technique [15]. For example, in the £ direction, 

I + h6*A 

= T i {l+h6\D\+k6{D- A )T^ 

= T((I - ahb Aj + h6\D + A )[I + ah\bAj\)~ l (I + ahD +• + h6{D A )T ^ 

= T i L A D A U A T[ l (4) 

where a is a coefficient appearing at the diagonal elements for the upwind differencing. 
In other words, a = 1 when first-order differencing is used for the ^-derivative, and 
3/2 when second-order differencing is used. Currently, a is set to be 7/6 since three- 
point first-order differencing is used for the ^-derivative to be consistent with fourth- 
order differencing in the right-hand side convective terms. The decomposition used in 
equation (4) is the approximate LDU factorization. This decomposition is more stable 
than simple LU factorization since the diagonal element always has the absolute value 
of the eigenvalues. 

In the solution process, an inversion in one direction consists of one scalar forward 
sweep and one scalar backward sweep. Thus, the LU-AD1 algorithm requires little 
additional memory and is easily vectorized. It may be noteworthy that LDU operator 
can be regarded as a symmetric Gauss-Seidel relaxation. 

The split eigenvalues in each direction, D A , D% and £>£, are modified corresponding 
to the implicit artificial dissipation terms and the viscous terms [5]. 
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Nonlinear Artificial Dissipation Model 

The fourth-order dissipation model has been widely used with central differencing in 
Navier-Stokes computations [3,5,13]. On the other hand, the second-order dissipation 
produced by proper upwinding works better than the fourth-order dissipation model in 
regions near discontinuities. Such dissipation, however, is not desirable for the rest of 
the flow field since it reduces the spatial accuracy to at most first order. Higher-order 
TVD upwind schemes can be constructed by introducing flux limiters [18,19]. 

A numerical flux for first-order upwind differencing of the convective terms can be 
written as 

F i+ 1/2 = (Fj + F i+l )/ 2 - dfj +l/i / 2 (5) 

where dFj +1 ^ 2 = \A J+ i/ 2 \{Qj+i - Qj)- Introducing a flux limiter $ to obtain second- 
order central differencing, equation (5) is rewritten as, 

F 1+in = (Fj + F Hl )l 2 -(I- i)dF? +w / 2- (6) 

To construct the present dissipation model, a simplified model of a flux is considered. 
The matrix, |A J+ x/ 2 | is replaced by the spectral radius = |C7| -I- cr^ so that the 
artificial dissipation terms can be simply evaluated; 

dFj + 1/2 = (^/J) y+1/2 (jg i+1 - JQ y ) (7) 

where the scaled spectral radius Oc/J at j + 1/2 is taken as an arithmetic average at 
j and j + 1, This simple form may badly diffuse solutions even with the limiter. A 
small constant e is introduced to the present model as well as conventional smoothing 
terms. A fourth-order dissipation term is also necessary to keep global stability to the 
central differencing. The final form of the present model is written as a combination of 
second- and fourth-order smoothing terms similar to an existing nonlinear model [17]. 

D t \zQ = ehVt((I - *, + i/2)d*7 + i/2 " *}dF? +l/2 ) (8) 

The limiter $ can be defined as a scalar function for each element: 


4>j+ 1/2 = mt'nmod(l,r + ,r ) 


(9) 
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where r + = dfj^/zldfj+x/^r- = dfj +3 / 2 /df } + 1/2 , and <f> and df are elements of $ 
and dF <T , respectively. 

The present dissipation model was implemented only in the £ direction where a strong 
discontinuity may appear. In the other directions, only the fourth-order dissipation 
term was used by setting $ to be I in equation (8). Only the fourth-order dissipation 
terms were used on the left-hand side of equation (3). 


RESULTS AND DISCUSSIONS 


Wing in Wind Tunnel 

The grid topology used here is of the C-H type with C-grid around the wing section 
and H-grid for the wing span. First, grid distributions in the y (spanwise) direction were 
determined. Then the geometry of each wing section was computed at each spanwise 
location by using the surface generator code E88F5 due to Sobieczky [12]. Finally, the 
C-grid system in the x-z plane was generated by using the two-dimensional algebraic 
technique for each wing section with outer boundaries fitted to the control-surface 
geometry. The grid was clustered near wall surfaces including the wing surface. The 
minimum spacing normal to surfaces was set to 5 x 10 -5 of the tunnel dimension which 
is 1 meter. The grid for each wing section at a constant spanwise location consists of 
201 x 51 points, and 75 planes are distributed in the spanwise direction for a total of 
768,825 points, of which 161 x 55 points are distributed on the wing surface. Figure 
1 shows three two-dimensional views of the grid system. Figure 2 shows the overall 
view of the wing and the tunnel walls with the typical pressure distributions. Two grid 
systems were generated, corresponding to the wing at 0 and 2° angle of attack. 

The test wing of this workshop corresponds to a practical design configuration which 
causes problems with grid generation. First, both the wing tip and trailing edge of 
this wind-tunnel test model have thickness. Therefore, the surface grids in the wing 
extension and wake cut were closed at the next grid point from the tip and trailing 
edge, respectively. Next, the test model has a large curvature near the wing tip and 
root. Thus, orthogonality of the grid system was enforced only near the wing surface 


7 


in the x-z plane, not near the wing tip and root-side wall (not in the x-y plane and y-z 
plane). 

The flow conditions consist of Mach number 0.82, Reynolds number 10 million based 
on the tunnel dimension lm, and 0 and 2° angles of attack. For turbulent flow cal- 
culations, the Baldwin-Lomax turbulence model modified for wing-wall configurations 
was used for the following four regions: the wing and the root-side wall, the wing ex- 
tension (flow-through condition) and the other side wall, the lower or upper wall and 
the root-side wall, and the lower or upper wall and the other side wall. Two types of 
turbulent-flow computations were performed. One assumed fully turbulent flow and 
the other assumed transition to turbulent flow following experimental data provided by 
the Workshop Organizer. In total, four cases were computed for the wing in a wind tun- 
nel. In the following computations, the subroutine for the turbulence model was called 
every 10 time steps to reduce CPU time, since the turbulence model computations 
represented 30% of the time required for one time step. 

At the inflow boundary, all variables were fixed to experimental data given by the 
inlet and exit conditions code E17F5 [12]. At the outflow boundary, the pressure was 
fixed by using the code E17F5 and the other variables were extrapolated. No-slip and 
adiabatic wall conditions were used on both the wing surface and wind-tunnel walls. 
The pressure on the surfaces was obtained from the normal momentum equation. The 
initial conditions were obtained by assuming the inflow profiles of flow quantities in the 
y-z plane to remain unchanged in the x direction. Therefore, boundary-layer profiles 
were assumed on the tunnel walls as the initial condition, while the slow-start technique 
was used for the wing surface. 

The CFL number was set to 20 initially and increased up to 400 during the first 1000 
iterations. The smoothing coefficient e started from 0.4 and decreased to 0.08 at the 
same time. To guarantee convergence, the computation was continued until the total 
time, nondimensionalized by the tunnel dimension and speed of sound, exceeded 4.5. 
This meant using over 6000 iterations because of the constraint of the time step at the 
wing root. Since CPU time/step/grid point was 10.0 n sec, the total CPU time was 
about 13 hr for one case. 

Figure 3 shows the computed results for the fully turbulent flow at 0° angle of attack. 
Figure 3a shows the coefficient of pressure at three spanwise locations, 2y/b = 0.2, 0.5 
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and 0.8. The pressure-contour plots on the wing upper surface and the corresponding 
surface-flow pattern are shown in Fig. 3b and c. Two types of separation can be 
observed. One is the trailing-edge separation and the other is the spiral-like separation 
near the wing root. The latter is due to the interaction between the boundary layers 
on the wing surface and the side wall. Figure 3d shows the pressure-contour plots on 
the lower and upper walls. Both plots coincide with each other because of 0° angle of 
attack. 

Figure 4 shows corresponding plots for fully turbulent flow at 2° angle of attack. In 
Fig. 4c, the surface-flow pattern indicates occurrence of shock-induced separation. This 
separation line is almost straight in the mid-span region, then turns to the trailing edge 
near the wing tip where the shock wave becomes weaker, and turns to the upstream. 
Flow separation at the wing root is also observed. The effect of angle of attack is 
apparent from the contour plots even on the lower and upper tunnel walls as shown in 
Fig. 4d. 

Two discrepancies can be pointed out in Figs. 3d and 4d. One is a steep rise of 
pressure near the inflow boundary, x = —0.3. The given inflow data specified only the 
boundary-layer velocity profile. Since density and temperature were assumed to be 
uniform at the inflow boundary, the numerical solutions changed rapidly to overcome 
this inconsistency. Density or temperature profiles will be required to obtain well- 
defined inflow data. (Pressure alone is not sufficient to define the boundary-layer 
properties.) 

The other discrepancy is the presence of kinks in the contours near the outflow 
boundary. There is a large pressure gradient in the y direction at the outflow boundary, 
x — 0.63, which is caused by the splitter plate (root-side wall) [20]. The pressure 
profiles obtained by experiment showed only two-dimensional profiles varying in the y 
direction and did not show the effect of changing angle of attack. But the numerical 
results indicate three dimensionality and the effect of changing angle of attack (for 
example at x = 0.5 in Figs. 3d and 4d), even though the effects are less than that 
caused by the splitter plate. More precise data will be required to specify the pressure 
profiles at the outflow boundary. 

Figures 5 shows corresponding plots for the transitional flow at 0° angle of attack. 
The pressure gradient on the wing surface is more moderate near the wing root when 
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compared with the fully turbulent flow case. Thus, no trailing-edge separation is ob- 
served in Fig. 5c. The pressure contours in Fig. 5b have kinks along the transition line 
where the turbulence model was turned on as a step function. 

Figure 6 shows the results of transitional flow at 2° angle of attack. Again, the 
pressure gradient on the wing surface is more moderate near the wing root in Fig. 6b 
and the flow is attached in Fig. 6c. However, the shock wave appears to be stronger 
and thus induces flow separation at mid span. 

Flow-separation patterns are different for the fully turbulent flow case and the tran- 
sitional flow case. Profiles of u-velocity are plotted in Fig. 7. The flow appears to 
accelerate in the laminar flow region for the transitional flow case, while the fully 
turbulent flow leads to a thick boundary layer. 

Isolated Wing in Free Air 

The grid system for the isolated wing was generated in the same way as for the wing 
in a wind tunnel. The spanwise grid distributions were taken from default of the code 
E88F5 as 40 points on the wing with ten points added for the wing extension and one 
additional point for implementing the symmetry condition at the wing root. Then C- 
grids of 201 x 51 points were generated at each spanwise station, with outer boundaries 
located 10 chord away. The minimum spacing normal to the wing surface was set to 
5 x 10 -5 of each chord length. The number of grid points on the wing surface was 
151 x 41 points. The total number of grid points was 201 x 51 x 51; that is, 522,801 
points. The same grid system was used for different angle-of-attack cases. 

The flow conditions consist of Mach number 0.82, Reynolds number 10 million based 
on the reference length of lm, and 0, 2, 5, and 8° angles of attack. The flow was assumed 
as fully turbulent. For the boundary conditions, free-stream values were specified at the 
far-field boundary. At the outflow boundary, the pressure was fixed to the free-stream 
value and the other variables were extrapolated. The slow-start technique was used to 
start iterations. 

Figure 8 shows the chordwise pressure distributions at 2 y/b = 0.2, 0.5 and 0.8 for 
four angle-of-attack cases, respectively. Figure 9 shows the pressure-contour plots on 
the upper surface of the wing. The suction peaks appear to be lower in the free-air 
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case compared with the corresponding wind-tunnel case. For the wind-tunnel cases, 
however, the suction peak becomes much higher on the lower surface (for example, 
compare Fig. 4a with Fig. 8b) since this is the pressure side in the wind tunnel. There- 
fore, the lift coefficient is larger for the free-air case even though the suction peaks 
are lower. Figure 10 shows the surface-flow patterns. No separation is observed at 0° 
angle-of-attack. The root separation observed in the wind-tunnel cases (because of the 
interaction of the side wall) did not occur for the free-air cases. At 2° angle-of-attack, 
the shock wave is weaker for the free-air case and thus the shock-induced separation is 
not observed. Large-scale separation is observed at higher angles of attack. 

The lift coefficient Cl as a function of angle of attack a is shown in Fig. 11. The curve 
shows the nonlinearity at higher angles of attack. The lift coefficient obtained for fully 
turbulent flow computations at 2° angle of attack in the wind tunnel is smaller than 
that obtained in free air because of the existence of the tunnel walls. The transitional 
turbulence model causes differences in lift coefficients of about 15% for the wind-tunnel 
simulations. This indicates that improvements of turbulence modeling are needed for 
aerodynamic predictions even at moderate angles of attack. Figure 12 shows the drag 
polar. The free-air case produces the highest drag at 0° angle of attack here, since the 
boundary layer on the side wall for the wind-tunnel cases was not taken into account 
in the present force calculation. The fully turbulent flow cases result in higher total 
drag than the transitional flow cases. 

CONCLUDING REMARKS 


The Navier-Stokes simulation of transonic flow past a test wing in a wind tunnel and 
in free air was carried out using the LU-ADI factorization algorithm. The Baldwin- 
Lomax model was used under the assumption of the fully turbulent flow. The effects 
of the tunnel walls are well demonstrated by the present simulations. Lift coefficients 
predicted by the wind-tunnel simulation were lower than those by the free-air simula- 
tion. 

The specified transition line from laminar to turbulent flow was also tried for the 
wind-tunnel case. This led to different results for the pressure distributions and the 
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separation patterns. In future, more work is required to simulate transitional flow 
correctly. 

In addition, the grid-generation technique was not robust enough to treat complex 
geometry accurately, such as large curvature around the wing tip and root, and finite 
thickness at the wing tip and trailing edge. Further development is needed in this area. 

Numerical results indicate that the accuracy of experimental data was insufficient to 
establish boundary conditions for computations. Both the inflow and outflow boundary 
profiles were insufficient to obtain smooth flow field profiles. More precise experimental 
data will be required for CFD validation. 
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Fig. 1 Grid distributions for wind-tunnel simulation. 
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Fig. 2 Overall view of geometry and pressure distributions on surfaces. 
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a) Chordwise pressure distributions. c) Surface-flow pattern on wing upper surface. 


Fig. 3 Computed results for a = 0° and fully turbulent flow case iAf,,,, = 0.82, Re = 1.0 X 10 7 . 











a) Chordwise pressure distributions. c) Surface-flow pattern on wing upper surface. 

Fig. 4 Computed results for a = 2° and fully turbulent flow case = 0.82, Re= 1 .0 X 1 0 7 . 
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d) Pressure contour plots on lower and upper walls. 
Fig. 4 Concluded. 
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Fig. 12 Lift coefficient versus drag coefficient; M 0 0 = 0.82, Re — 1.0 x 10 7 . 
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